<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of fast_proj_b2</title>
  <meta name="keywords" content="fast_proj_b2">
  <meta name="description" content="PROJ_B2 - Projection onto a L2-ball">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../../index.html">Home</a> &gt;  <a href="../../index.html">CO4OI</a> &gt; <a href="#">Solvers</a> &gt; <a href="index.html">NM_solvers</a> &gt; fast_proj_b2.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../../index.html"><img alt="<" border="0" src="../../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for CO4OI/Solvers/NM_solvers&nbsp;<img alt=">" border="0" src="../../../right.png"></a></td></tr></table>-->

<h1>fast_proj_b2
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>PROJ_B2 - Projection onto a L2-ball</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>function [sol, u] = fast_proj_b2(x, param) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> PROJ_B2 - Projection onto a L2-ball

 sol = fast_proj_b2(x, param) solves:

   min_{z} ||x - z||_2^2   s.t.  ||y - A z||_2 &lt; epsilon

 param is a Matlab structure containing the following fields:

   - y: measurements (default: 0).

   - A: Forward operator (default: Id).

   - At: Adjoint operator (default: Id).

   - epsilon: Radius of the L2 ball (default = 1e-3).

   - tight: 1 if A is a tight frame or 0 if not (default = 1)

   - nu: bound on the norm of the operator A, i.e.
       ||A x||^2 &lt;= nu * ||x||^2 (default: 1)

   - tol: tolerance for the projection onto the L2 ball. The algorithms
   stops if
       epsilon/(1-tol) &lt;= ||y - A z||_2 &lt;= epsilon/(1+tol)
   (default: 1e-3)

   - max_iter: max. nb. of iterations (default: 200).

   - verbose: 0 no log, 1 a summary at convergence, 2 print main
   steps (default: 1)



 References:
 [1] M.J. Fadili and J-L. Starck, &quot;Monotone operator splitting for
 optimization problems in sparse recovery&quot; , IEEE ICIP, Cairo,
 Egypt, 2009.
 [2] Amir Beck and Marc Teboulle, &quot;A Fast Iterative Shrinkage-Thresholding
 Algorithm for Linear Inverse Problems&quot;,  SIAM Journal on Imaging Sciences
 2 (2009), no. 1, 183--202.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../../matlabicon.gif)">
</ul>
This function is called by:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="solve_nm.html" class="code" title="function sol = solve_nm(y, param)">solve_nm</a>	</li></ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [sol, u] = fast_proj_b2(x, param)</a>
0002 <span class="comment">% PROJ_B2 - Projection onto a L2-ball</span>
0003 <span class="comment">%</span>
0004 <span class="comment">% sol = fast_proj_b2(x, param) solves:</span>
0005 <span class="comment">%</span>
0006 <span class="comment">%   min_{z} ||x - z||_2^2   s.t.  ||y - A z||_2 &lt; epsilon</span>
0007 <span class="comment">%</span>
0008 <span class="comment">% param is a Matlab structure containing the following fields:</span>
0009 <span class="comment">%</span>
0010 <span class="comment">%   - y: measurements (default: 0).</span>
0011 <span class="comment">%</span>
0012 <span class="comment">%   - A: Forward operator (default: Id).</span>
0013 <span class="comment">%</span>
0014 <span class="comment">%   - At: Adjoint operator (default: Id).</span>
0015 <span class="comment">%</span>
0016 <span class="comment">%   - epsilon: Radius of the L2 ball (default = 1e-3).</span>
0017 <span class="comment">%</span>
0018 <span class="comment">%   - tight: 1 if A is a tight frame or 0 if not (default = 1)</span>
0019 <span class="comment">%</span>
0020 <span class="comment">%   - nu: bound on the norm of the operator A, i.e.</span>
0021 <span class="comment">%       ||A x||^2 &lt;= nu * ||x||^2 (default: 1)</span>
0022 <span class="comment">%</span>
0023 <span class="comment">%   - tol: tolerance for the projection onto the L2 ball. The algorithms</span>
0024 <span class="comment">%   stops if</span>
0025 <span class="comment">%       epsilon/(1-tol) &lt;= ||y - A z||_2 &lt;= epsilon/(1+tol)</span>
0026 <span class="comment">%   (default: 1e-3)</span>
0027 <span class="comment">%</span>
0028 <span class="comment">%   - max_iter: max. nb. of iterations (default: 200).</span>
0029 <span class="comment">%</span>
0030 <span class="comment">%   - verbose: 0 no log, 1 a summary at convergence, 2 print main</span>
0031 <span class="comment">%   steps (default: 1)</span>
0032 <span class="comment">%</span>
0033 <span class="comment">%</span>
0034 <span class="comment">%</span>
0035 <span class="comment">% References:</span>
0036 <span class="comment">% [1] M.J. Fadili and J-L. Starck, &quot;Monotone operator splitting for</span>
0037 <span class="comment">% optimization problems in sparse recovery&quot; , IEEE ICIP, Cairo,</span>
0038 <span class="comment">% Egypt, 2009.</span>
0039 <span class="comment">% [2] Amir Beck and Marc Teboulle, &quot;A Fast Iterative Shrinkage-Thresholding</span>
0040 <span class="comment">% Algorithm for Linear Inverse Problems&quot;,  SIAM Journal on Imaging Sciences</span>
0041 <span class="comment">% 2 (2009), no. 1, 183--202.</span>
0042 <span class="comment">%</span>
0043 
0044 
0045 <span class="comment">% Optional input arguments</span>
0046 <span class="keyword">if</span> ~isfield(param, <span class="string">'y'</span>), param.y = 0; <span class="keyword">end</span>
0047 <span class="keyword">if</span> ~isfield(param, <span class="string">'A'</span>), param.A = @(x) x; param.At = @(x) x; <span class="keyword">end</span>
0048 <span class="keyword">if</span> ~isfield(param, <span class="string">'epsilon'</span>), param.epsilon = 1e-3; <span class="keyword">end</span>
0049 <span class="keyword">if</span> ~isfield(param, <span class="string">'tight'</span>), param.tight = 1; <span class="keyword">end</span>
0050 <span class="keyword">if</span> ~isfield(param, <span class="string">'tol'</span>), param.tol = 1e-3; <span class="keyword">end</span>
0051 <span class="keyword">if</span> ~isfield(param, <span class="string">'verbose'</span>), param.verbose = 1; <span class="keyword">end</span>
0052 <span class="keyword">if</span> ~isfield(param, <span class="string">'nu'</span>), param.nu = 1; <span class="keyword">end</span>
0053 <span class="keyword">if</span> ~isfield(param, <span class="string">'max_iter'</span>), param.max_iter = 200; <span class="keyword">end</span>
0054 <span class="keyword">if</span> ~isfield(param, <span class="string">'u'</span>), param.u = zeros(size(param.y)); <span class="keyword">end</span>
0055 <span class="keyword">if</span> ~isfield(param, <span class="string">'pos'</span>), param.pos = 0; <span class="keyword">end</span>
0056 <span class="keyword">if</span> ~isfield(param, <span class="string">'real'</span>), param.real = 0; <span class="keyword">end</span>
0057 
0058 <span class="comment">% Useful functions for the projection</span>
0059 sc = @(z) z*min(param.epsilon/norm(z(:)), 1); <span class="comment">% scaling</span>
0060 
0061 <span class="comment">% Projection</span>
0062 
0063 <span class="comment">% TIGHT FRAME CASE</span>
0064 <span class="comment">% In NM algorithm this is the used case, the general (non-tight) case is</span>
0065 <span class="comment">% included for completeness</span>
0066 <span class="keyword">if</span> (param.tight &amp;&amp; ~(param.pos||param.real))
0067     
0068     temp = param.A(x) - param.y;
0069     sol = x + 1/param.nu * param.At(sc(temp)-temp);
0070     crit_B2 = <span class="string">'TOL_EPS'</span>; iter = 0;
0071     u = 0;
0072     
0073 <span class="keyword">else</span> <span class="comment">% NON TIGHT FRAME CASE</span>
0074     
0075     <span class="comment">% Initializations</span>
0076     sol = x; u = param.u; v = u;
0077     iter = 1; true = 1; told = 1;
0078     
0079     <span class="comment">% Tolerance onto the L2 ball</span>
0080     epsilon_low = param.epsilon_low;
0081     epsilon_up = param.epsilon_up;
0082     
0083     <span class="comment">% Check if we are in the L2 ball</span>
0084     dummy = param.A(sol);
0085     norm_res = norm(param.y(:)-dummy(:), 2);
0086     <span class="keyword">if</span> norm_res &lt;= epsilon_up
0087         crit_B2 = <span class="string">'TOL_EPS'</span>; true = 0;
0088     <span class="keyword">end</span>
0089     
0090     <span class="comment">% Projection onto the L2-ball</span>
0091     <span class="comment">% Init</span>
0092     <span class="keyword">if</span> param.verbose &gt; 1
0093         fprintf(<span class="string">'  Proj. B2:\n'</span>);
0094     <span class="keyword">end</span>
0095     <span class="keyword">while</span> true
0096         
0097         <span class="comment">% Residual</span>
0098         res = param.A(sol) - param.y; norm_res = norm(res(:), 2);
0099         
0100         <span class="comment">% Scaling for the projection</span>
0101         res = u*param.nu + res; norm_proj = norm(res(:), 2);
0102         
0103         <span class="comment">% Log</span>
0104         <span class="keyword">if</span> param.verbose&gt;1
0105             fprintf(<span class="string">'   Iter %i, epsilon = %e, ||y - Ax||_2 = %e\n'</span>, <span class="keyword">...</span>
0106                 iter, param.epsilon, norm_res);
0107         <span class="keyword">end</span>
0108         
0109         <span class="comment">% Stopping criterion</span>
0110         <span class="keyword">if</span> (norm_res&gt;=epsilon_low &amp;&amp; norm_res&lt;=epsilon_up)
0111             crit_B2 = <span class="string">'TOL_EPS'</span>; <span class="keyword">break</span>;
0112         <span class="keyword">elseif</span> iter &gt;= param.max_iter
0113             crit_B2 = <span class="string">'MAX_IT'</span>; <span class="keyword">break</span>;
0114         <span class="keyword">end</span>
0115         
0116         <span class="comment">% Projection onto the L2 ball and FISTA update</span>
0117         t = (1+sqrt(1+4*told^2))/2;
0118         ratio = min(1, param.epsilon/norm_proj);
0119         u = v;
0120         v = 1/param.nu * (res - res*ratio);
0121         u = v + (told-1)/t * (v - u);
0122         
0123         <span class="comment">% Current estimate</span>
0124         sol = x - param.At(u);
0125         
0126         <span class="comment">%Projection onto the non-negative orthant (positivity constraint)</span>
0127         
0128         <span class="keyword">if</span> (param.pos)
0129             sol=real(sol);
0130             sol(sol&lt;0)=0;
0131             
0132         <span class="keyword">end</span>
0133         
0134         <span class="comment">%Projection onto the real orthant (reality constraint)</span>
0135         
0136         <span class="keyword">if</span> (param.real)
0137             sol=real(sol);
0138             
0139         <span class="keyword">end</span>
0140         
0141         
0142         
0143         <span class="comment">% Update number of iteration</span>
0144         told = t;
0145         
0146         <span class="comment">% Update number of iterations</span>
0147         iter = iter + 1;
0148         
0149     <span class="keyword">end</span>
0150 <span class="keyword">end</span>
0151 
0152 <span class="comment">% Log after the projection onto the L2-ball</span>
0153 <span class="keyword">if</span> param.verbose &gt;= 1
0154     temp = param.A(sol);
0155     fprintf([<span class="string">'  Proj. B2: epsilon = %e, ||y-Ax||_2 = %e,'</span>, <span class="keyword">...</span>
0156         <span class="string">' %s, iter = %i\n'</span>], param.epsilon, norm(param.y(:)-temp(:)), <span class="keyword">...</span>
0157         crit_B2, iter);
0158 <span class="keyword">end</span>
0159 
0160 
0161 <span class="keyword">end</span></pre></div>
<hr><address>Generated on Thu 18-Jul-2013 19:25:15 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>